In [1]:
using XSim
using Distributions
using StatsFuns
using Gadfly
using Random
using Printf
using StatsBase
using LinearAlgebra
using JWAS
using DataFrames
using DelimitedFiles
using Cairo
using Fontconfig
using GLM
In [2]:
include("Functions.jl")
Out[2]:
SNP_SNP_GWAS2 (generic function with 2 methods)
In [3]:
Random.seed!(95616)
numChr,numLoci,chrLength,mutRate = 1,3020,1.0,2.5e-5  # Mutation rate was 2.5*10−5 during the historical generations
locusInt  = chrLength/numLoci
mapPos   = collect(0:locusInt:(chrLength-0.0001))  #Uniform position
geneFreq   = fill(0.5,numLoci)  # initial gene frequency 0.5
XSim.build_genome(numChr,chrLength,numLoci,geneFreq,mapPos,mutRate)
In [4]:
sireSizeFounder = 200
damSizeFounder = 200
sires = sampleFounders(sireSizeFounder)
dams  = sampleFounders(damSizeFounder);
Sampling 200 animals into base population.
Sampling 200 animals into base population.
In [5]:
ngen,popSize = 500,400
sires1,dams1,gen1 = sampleRan(popSize, ngen, sires, dams);
Generation     2: sampling   200 males and   200 females
Generation     3: sampling   200 males and   200 females
Generation     4: sampling   200 males and   200 females
Generation     5: sampling   200 males and   200 females
Generation     6: sampling   200 males and   200 females
Generation     7: sampling   200 males and   200 females
Generation     8: sampling   200 males and   200 females
Generation     9: sampling   200 males and   200 females
Generation    10: sampling   200 males and   200 females
Generation    11: sampling   200 males and   200 females
Generation    12: sampling   200 males and   200 females
Generation    13: sampling   200 males and   200 females
Generation    14: sampling   200 males and   200 females
Generation    15: sampling   200 males and   200 females
Generation    16: sampling   200 males and   200 females
Generation    17: sampling   200 males and   200 females
Generation    18: sampling   200 males and   200 females
Generation    19: sampling   200 males and   200 females
Generation    20: sampling   200 males and   200 females
Generation    21: sampling   200 males and   200 females
Generation    22: sampling   200 males and   200 females
Generation    23: sampling   200 males and   200 females
Generation    24: sampling   200 males and   200 females
Generation    25: sampling   200 males and   200 females
Generation    26: sampling   200 males and   200 females
Generation    27: sampling   200 males and   200 females
Generation    28: sampling   200 males and   200 females
Generation    29: sampling   200 males and   200 females
Generation    30: sampling   200 males and   200 females
Generation    31: sampling   200 males and   200 females
Generation    32: sampling   200 males and   200 females
Generation    33: sampling   200 males and   200 females
Generation    34: sampling   200 males and   200 females
Generation    35: sampling   200 males and   200 females
Generation    36: sampling   200 males and   200 females
Generation    37: sampling   200 males and   200 females
Generation    38: sampling   200 males and   200 females
Generation    39: sampling   200 males and   200 females
Generation    40: sampling   200 males and   200 females
Generation    41: sampling   200 males and   200 females
Generation    42: sampling   200 males and   200 females
Generation    43: sampling   200 males and   200 females
Generation    44: sampling   200 males and   200 females
Generation    45: sampling   200 males and   200 females
Generation    46: sampling   200 males and   200 females
Generation    47: sampling   200 males and   200 females
Generation    48: sampling   200 males and   200 females
Generation    49: sampling   200 males and   200 females
Generation    50: sampling   200 males and   200 females
Generation    51: sampling   200 males and   200 females
Generation    52: sampling   200 males and   200 females
Generation    53: sampling   200 males and   200 females
Generation    54: sampling   200 males and   200 females
Generation    55: sampling   200 males and   200 females
Generation    56: sampling   200 males and   200 females
Generation    57: sampling   200 males and   200 females
Generation    58: sampling   200 males and   200 females
Generation    59: sampling   200 males and   200 females
Generation    60: sampling   200 males and   200 females
Generation    61: sampling   200 males and   200 females
Generation    62: sampling   200 males and   200 females
Generation    63: sampling   200 males and   200 females
Generation    64: sampling   200 males and   200 females
Generation    65: sampling   200 males and   200 females
Generation    66: sampling   200 males and   200 females
Generation    67: sampling   200 males and   200 females
Generation    68: sampling   200 males and   200 females
Generation    69: sampling   200 males and   200 females
Generation    70: sampling   200 males and   200 females
Generation    71: sampling   200 males and   200 females
Generation    72: sampling   200 males and   200 females
Generation    73: sampling   200 males and   200 females
Generation    74: sampling   200 males and   200 females
Generation    75: sampling   200 males and   200 females
Generation    76: sampling   200 males and   200 females
Generation    77: sampling   200 males and   200 females
Generation    78: sampling   200 males and   200 females
Generation    79: sampling   200 males and   200 females
Generation    80: sampling   200 males and   200 females
Generation    81: sampling   200 males and   200 females
Generation    82: sampling   200 males and   200 females
Generation    83: sampling   200 males and   200 females
Generation    84: sampling   200 males and   200 females
Generation    85: sampling   200 males and   200 females
Generation    86: sampling   200 males and   200 females
Generation    87: sampling   200 males and   200 females
Generation    88: sampling   200 males and   200 females
Generation    89: sampling   200 males and   200 females
Generation    90: sampling   200 males and   200 females
Generation    91: sampling   200 males and   200 females
Generation    92: sampling   200 males and   200 females
Generation    93: sampling   200 males and   200 females
Generation    94: sampling   200 males and   200 females
Generation    95: sampling   200 males and   200 females
Generation    96: sampling   200 males and   200 females
Generation    97: sampling   200 males and   200 females
Generation    98: sampling   200 males and   200 females
Generation    99: sampling   200 males and   200 females
Generation   100: sampling   200 males and   200 females
Generation   101: sampling   200 males and   200 females
Generation   102: sampling   200 males and   200 females
Generation   103: sampling   200 males and   200 females
Generation   104: sampling   200 males and   200 females
Generation   105: sampling   200 males and   200 females
Generation   106: sampling   200 males and   200 females
Generation   107: sampling   200 males and   200 females
Generation   108: sampling   200 males and   200 females
Generation   109: sampling   200 males and   200 females
Generation   110: sampling   200 males and   200 females
Generation   111: sampling   200 males and   200 females
Generation   112: sampling   200 males and   200 females
Generation   113: sampling   200 males and   200 females
Generation   114: sampling   200 males and   200 females
Generation   115: sampling   200 males and   200 females
Generation   116: sampling   200 males and   200 females
Generation   117: sampling   200 males and   200 females
Generation   118: sampling   200 males and   200 females
Generation   119: sampling   200 males and   200 females
Generation   120: sampling   200 males and   200 females
Generation   121: sampling   200 males and   200 females
Generation   122: sampling   200 males and   200 females
Generation   123: sampling   200 males and   200 females
Generation   124: sampling   200 males and   200 females
Generation   125: sampling   200 males and   200 females
Generation   126: sampling   200 males and   200 females
Generation   127: sampling   200 males and   200 females
Generation   128: sampling   200 males and   200 females
Generation   129: sampling   200 males and   200 females
Generation   130: sampling   200 males and   200 females
Generation   131: sampling   200 males and   200 females
Generation   132: sampling   200 males and   200 females
Generation   133: sampling   200 males and   200 females
Generation   134: sampling   200 males and   200 females
Generation   135: sampling   200 males and   200 females
Generation   136: sampling   200 males and   200 females
Generation   137: sampling   200 males and   200 females
Generation   138: sampling   200 males and   200 females
Generation   139: sampling   200 males and   200 females
Generation   140: sampling   200 males and   200 females
Generation   141: sampling   200 males and   200 females
Generation   142: sampling   200 males and   200 females
Generation   143: sampling   200 males and   200 females
Generation   144: sampling   200 males and   200 females
Generation   145: sampling   200 males and   200 females
Generation   146: sampling   200 males and   200 females
Generation   147: sampling   200 males and   200 females
Generation   148: sampling   200 males and   200 females
Generation   149: sampling   200 males and   200 females
Generation   150: sampling   200 males and   200 females
Generation   151: sampling   200 males and   200 females
Generation   152: sampling   200 males and   200 females
Generation   153: sampling   200 males and   200 females
Generation   154: sampling   200 males and   200 females
Generation   155: sampling   200 males and   200 females
Generation   156: sampling   200 males and   200 females
Generation   157: sampling   200 males and   200 females
Generation   158: sampling   200 males and   200 females
Generation   159: sampling   200 males and   200 females
Generation   160: sampling   200 males and   200 females
Generation   161: sampling   200 males and   200 females
Generation   162: sampling   200 males and   200 females
Generation   163: sampling   200 males and   200 females
Generation   164: sampling   200 males and   200 females
Generation   165: sampling   200 males and   200 females
Generation   166: sampling   200 males and   200 females
Generation   167: sampling   200 males and   200 females
Generation   168: sampling   200 males and   200 females
Generation   169: sampling   200 males and   200 females
Generation   170: sampling   200 males and   200 females
Generation   171: sampling   200 males and   200 females
Generation   172: sampling   200 males and   200 females
Generation   173: sampling   200 males and   200 females
Generation   174: sampling   200 males and   200 females
Generation   175: sampling   200 males and   200 females
Generation   176: sampling   200 males and   200 females
Generation   177: sampling   200 males and   200 females
Generation   178: sampling   200 males and   200 females
Generation   179: sampling   200 males and   200 females
Generation   180: sampling   200 males and   200 females
Generation   181: sampling   200 males and   200 females
Generation   182: sampling   200 males and   200 females
Generation   183: sampling   200 males and   200 females
Generation   184: sampling   200 males and   200 females
Generation   185: sampling   200 males and   200 females
Generation   186: sampling   200 males and   200 females
Generation   187: sampling   200 males and   200 females
Generation   188: sampling   200 males and   200 females
Generation   189: sampling   200 males and   200 females
Generation   190: sampling   200 males and   200 females
Generation   191: sampling   200 males and   200 females
Generation   192: sampling   200 males and   200 females
Generation   193: sampling   200 males and   200 females
Generation   194: sampling   200 males and   200 females
Generation   195: sampling   200 males and   200 females
Generation   196: sampling   200 males and   200 females
Generation   197: sampling   200 males and   200 females
Generation   198: sampling   200 males and   200 females
Generation   199: sampling   200 males and   200 females
Generation   200: sampling   200 males and   200 females
Generation   201: sampling   200 males and   200 females
Generation   202: sampling   200 males and   200 females
Generation   203: sampling   200 males and   200 females
Generation   204: sampling   200 males and   200 females
Generation   205: sampling   200 males and   200 females
Generation   206: sampling   200 males and   200 females
Generation   207: sampling   200 males and   200 females
Generation   208: sampling   200 males and   200 females
Generation   209: sampling   200 males and   200 females
Generation   210: sampling   200 males and   200 females
Generation   211: sampling   200 males and   200 females
Generation   212: sampling   200 males and   200 females
Generation   213: sampling   200 males and   200 females
Generation   214: sampling   200 males and   200 females
Generation   215: sampling   200 males and   200 females
Generation   216: sampling   200 males and   200 females
Generation   217: sampling   200 males and   200 females
Generation   218: sampling   200 males and   200 females
Generation   219: sampling   200 males and   200 females
Generation   220: sampling   200 males and   200 females
Generation   221: sampling   200 males and   200 females
Generation   222: sampling   200 males and   200 females
Generation   223: sampling   200 males and   200 females
Generation   224: sampling   200 males and   200 females
Generation   225: sampling   200 males and   200 females
Generation   226: sampling   200 males and   200 females
Generation   227: sampling   200 males and   200 females
Generation   228: sampling   200 males and   200 females
Generation   229: sampling   200 males and   200 females
Generation   230: sampling   200 males and   200 females
Generation   231: sampling   200 males and   200 females
Generation   232: sampling   200 males and   200 females
Generation   233: sampling   200 males and   200 females
Generation   234: sampling   200 males and   200 females
Generation   235: sampling   200 males and   200 females
Generation   236: sampling   200 males and   200 females
Generation   237: sampling   200 males and   200 females
Generation   238: sampling   200 males and   200 females
Generation   239: sampling   200 males and   200 females
Generation   240: sampling   200 males and   200 females
Generation   241: sampling   200 males and   200 females
Generation   242: sampling   200 males and   200 females
Generation   243: sampling   200 males and   200 females
Generation   244: sampling   200 males and   200 females
Generation   245: sampling   200 males and   200 females
Generation   246: sampling   200 males and   200 females
Generation   247: sampling   200 males and   200 females
Generation   248: sampling   200 males and   200 females
Generation   249: sampling   200 males and   200 females
Generation   250: sampling   200 males and   200 females
Generation   251: sampling   200 males and   200 females
Generation   252: sampling   200 males and   200 females
Generation   253: sampling   200 males and   200 females
Generation   254: sampling   200 males and   200 females
Generation   255: sampling   200 males and   200 females
Generation   256: sampling   200 males and   200 females
Generation   257: sampling   200 males and   200 females
Generation   258: sampling   200 males and   200 females
Generation   259: sampling   200 males and   200 females
Generation   260: sampling   200 males and   200 females
Generation   261: sampling   200 males and   200 females
Generation   262: sampling   200 males and   200 females
Generation   263: sampling   200 males and   200 females
Generation   264: sampling   200 males and   200 females
Generation   265: sampling   200 males and   200 females
Generation   266: sampling   200 males and   200 females
Generation   267: sampling   200 males and   200 females
Generation   268: sampling   200 males and   200 females
Generation   269: sampling   200 males and   200 females
Generation   270: sampling   200 males and   200 females
Generation   271: sampling   200 males and   200 females
Generation   272: sampling   200 males and   200 females
Generation   273: sampling   200 males and   200 females
Generation   274: sampling   200 males and   200 females
Generation   275: sampling   200 males and   200 females
Generation   276: sampling   200 males and   200 females
Generation   277: sampling   200 males and   200 females
Generation   278: sampling   200 males and   200 females
Generation   279: sampling   200 males and   200 females
Generation   280: sampling   200 males and   200 females
Generation   281: sampling   200 males and   200 females
Generation   282: sampling   200 males and   200 females
Generation   283: sampling   200 males and   200 females
Generation   284: sampling   200 males and   200 females
Generation   285: sampling   200 males and   200 females
Generation   286: sampling   200 males and   200 females
Generation   287: sampling   200 males and   200 females
Generation   288: sampling   200 males and   200 females
Generation   289: sampling   200 males and   200 females
Generation   290: sampling   200 males and   200 females
Generation   291: sampling   200 males and   200 females
Generation   292: sampling   200 males and   200 females
Generation   293: sampling   200 males and   200 females
Generation   294: sampling   200 males and   200 females
Generation   295: sampling   200 males and   200 females
Generation   296: sampling   200 males and   200 females
Generation   297: sampling   200 males and   200 females
Generation   298: sampling   200 males and   200 females
Generation   299: sampling   200 males and   200 females
Generation   300: sampling   200 males and   200 females
Generation   301: sampling   200 males and   200 females
Generation   302: sampling   200 males and   200 females
Generation   303: sampling   200 males and   200 females
Generation   304: sampling   200 males and   200 females
Generation   305: sampling   200 males and   200 females
Generation   306: sampling   200 males and   200 females
Generation   307: sampling   200 males and   200 females
Generation   308: sampling   200 males and   200 females
Generation   309: sampling   200 males and   200 females
Generation   310: sampling   200 males and   200 females
Generation   311: sampling   200 males and   200 females
Generation   312: sampling   200 males and   200 females
Generation   313: sampling   200 males and   200 females
Generation   314: sampling   200 males and   200 females
Generation   315: sampling   200 males and   200 females
Generation   316: sampling   200 males and   200 females
Generation   317: sampling   200 males and   200 females
Generation   318: sampling   200 males and   200 females
Generation   319: sampling   200 males and   200 females
Generation   320: sampling   200 males and   200 females
Generation   321: sampling   200 males and   200 females
Generation   322: sampling   200 males and   200 females
Generation   323: sampling   200 males and   200 females
Generation   324: sampling   200 males and   200 females
Generation   325: sampling   200 males and   200 females
Generation   326: sampling   200 males and   200 females
Generation   327: sampling   200 males and   200 females
Generation   328: sampling   200 males and   200 females
Generation   329: sampling   200 males and   200 females
Generation   330: sampling   200 males and   200 females
Generation   331: sampling   200 males and   200 females
Generation   332: sampling   200 males and   200 females
Generation   333: sampling   200 males and   200 females
Generation   334: sampling   200 males and   200 females
Generation   335: sampling   200 males and   200 females
Generation   336: sampling   200 males and   200 females
Generation   337: sampling   200 males and   200 females
Generation   338: sampling   200 males and   200 females
Generation   339: sampling   200 males and   200 females
Generation   340: sampling   200 males and   200 females
Generation   341: sampling   200 males and   200 females
Generation   342: sampling   200 males and   200 females
Generation   343: sampling   200 males and   200 females
Generation   344: sampling   200 males and   200 females
Generation   345: sampling   200 males and   200 females
Generation   346: sampling   200 males and   200 females
Generation   347: sampling   200 males and   200 females
Generation   348: sampling   200 males and   200 females
Generation   349: sampling   200 males and   200 females
Generation   350: sampling   200 males and   200 females
Generation   351: sampling   200 males and   200 females
Generation   352: sampling   200 males and   200 females
Generation   353: sampling   200 males and   200 females
Generation   354: sampling   200 males and   200 females
Generation   355: sampling   200 males and   200 females
Generation   356: sampling   200 males and   200 females
Generation   357: sampling   200 males and   200 females
Generation   358: sampling   200 males and   200 females
Generation   359: sampling   200 males and   200 females
Generation   360: sampling   200 males and   200 females
Generation   361: sampling   200 males and   200 females
Generation   362: sampling   200 males and   200 females
Generation   363: sampling   200 males and   200 females
Generation   364: sampling   200 males and   200 females
Generation   365: sampling   200 males and   200 females
Generation   366: sampling   200 males and   200 females
Generation   367: sampling   200 males and   200 females
Generation   368: sampling   200 males and   200 females
Generation   369: sampling   200 males and   200 females
Generation   370: sampling   200 males and   200 females
Generation   371: sampling   200 males and   200 females
Generation   372: sampling   200 males and   200 females
Generation   373: sampling   200 males and   200 females
Generation   374: sampling   200 males and   200 females
Generation   375: sampling   200 males and   200 females
Generation   376: sampling   200 males and   200 females
Generation   377: sampling   200 males and   200 females
Generation   378: sampling   200 males and   200 females
Generation   379: sampling   200 males and   200 females
Generation   380: sampling   200 males and   200 females
Generation   381: sampling   200 males and   200 females
Generation   382: sampling   200 males and   200 females
Generation   383: sampling   200 males and   200 females
Generation   384: sampling   200 males and   200 females
Generation   385: sampling   200 males and   200 females
Generation   386: sampling   200 males and   200 females
Generation   387: sampling   200 males and   200 females
Generation   388: sampling   200 males and   200 females
Generation   389: sampling   200 males and   200 females
Generation   390: sampling   200 males and   200 females
Generation   391: sampling   200 males and   200 females
Generation   392: sampling   200 males and   200 females
Generation   393: sampling   200 males and   200 females
Generation   394: sampling   200 males and   200 females
Generation   395: sampling   200 males and   200 females
Generation   396: sampling   200 males and   200 females
Generation   397: sampling   200 males and   200 females
Generation   398: sampling   200 males and   200 females
Generation   399: sampling   200 males and   200 females
Generation   400: sampling   200 males and   200 females
Generation   401: sampling   200 males and   200 females
Generation   402: sampling   200 males and   200 females
Generation   403: sampling   200 males and   200 females
Generation   404: sampling   200 males and   200 females
Generation   405: sampling   200 males and   200 females
Generation   406: sampling   200 males and   200 females
Generation   407: sampling   200 males and   200 females
Generation   408: sampling   200 males and   200 females
Generation   409: sampling   200 males and   200 females
Generation   410: sampling   200 males and   200 females
Generation   411: sampling   200 males and   200 females
Generation   412: sampling   200 males and   200 females
Generation   413: sampling   200 males and   200 females
Generation   414: sampling   200 males and   200 females
Generation   415: sampling   200 males and   200 females
Generation   416: sampling   200 males and   200 females
Generation   417: sampling   200 males and   200 females
Generation   418: sampling   200 males and   200 females
Generation   419: sampling   200 males and   200 females
Generation   420: sampling   200 males and   200 females
Generation   421: sampling   200 males and   200 females
Generation   422: sampling   200 males and   200 females
Generation   423: sampling   200 males and   200 females
Generation   424: sampling   200 males and   200 females
Generation   425: sampling   200 males and   200 females
Generation   426: sampling   200 males and   200 females
Generation   427: sampling   200 males and   200 females
Generation   428: sampling   200 males and   200 females
Generation   429: sampling   200 males and   200 females
Generation   430: sampling   200 males and   200 females
Generation   431: sampling   200 males and   200 females
Generation   432: sampling   200 males and   200 females
Generation   433: sampling   200 males and   200 females
Generation   434: sampling   200 males and   200 females
Generation   435: sampling   200 males and   200 females
Generation   436: sampling   200 males and   200 females
Generation   437: sampling   200 males and   200 females
Generation   438: sampling   200 males and   200 females
Generation   439: sampling   200 males and   200 females
Generation   440: sampling   200 males and   200 females
Generation   441: sampling   200 males and   200 females
Generation   442: sampling   200 males and   200 females
Generation   443: sampling   200 males and   200 females
Generation   444: sampling   200 males and   200 females
Generation   445: sampling   200 males and   200 females
Generation   446: sampling   200 males and   200 females
Generation   447: sampling   200 males and   200 females
Generation   448: sampling   200 males and   200 females
Generation   449: sampling   200 males and   200 females
Generation   450: sampling   200 males and   200 females
Generation   451: sampling   200 males and   200 females
Generation   452: sampling   200 males and   200 females
Generation   453: sampling   200 males and   200 females
Generation   454: sampling   200 males and   200 females
Generation   455: sampling   200 males and   200 females
Generation   456: sampling   200 males and   200 females
Generation   457: sampling   200 males and   200 females
Generation   458: sampling   200 males and   200 females
Generation   459: sampling   200 males and   200 females
Generation   460: sampling   200 males and   200 females
Generation   461: sampling   200 males and   200 females
Generation   462: sampling   200 males and   200 females
Generation   463: sampling   200 males and   200 females
Generation   464: sampling   200 males and   200 females
Generation   465: sampling   200 males and   200 females
Generation   466: sampling   200 males and   200 females
Generation   467: sampling   200 males and   200 females
Generation   468: sampling   200 males and   200 females
Generation   469: sampling   200 males and   200 females
Generation   470: sampling   200 males and   200 females
Generation   471: sampling   200 males and   200 females
Generation   472: sampling   200 males and   200 females
Generation   473: sampling   200 males and   200 females
Generation   474: sampling   200 males and   200 females
Generation   475: sampling   200 males and   200 females
Generation   476: sampling   200 males and   200 females
Generation   477: sampling   200 males and   200 females
Generation   478: sampling   200 males and   200 females
Generation   479: sampling   200 males and   200 females
Generation   480: sampling   200 males and   200 females
Generation   481: sampling   200 males and   200 females
Generation   482: sampling   200 males and   200 females
Generation   483: sampling   200 males and   200 females
Generation   484: sampling   200 males and   200 females
Generation   485: sampling   200 males and   200 females
Generation   486: sampling   200 males and   200 females
Generation   487: sampling   200 males and   200 females
Generation   488: sampling   200 males and   200 females
Generation   489: sampling   200 males and   200 females
Generation   490: sampling   200 males and   200 females
Generation   491: sampling   200 males and   200 females
Generation   492: sampling   200 males and   200 females
Generation   493: sampling   200 males and   200 females
Generation   494: sampling   200 males and   200 females
Generation   495: sampling   200 males and   200 females
Generation   496: sampling   200 males and   200 females
Generation   497: sampling   200 males and   200 females
Generation   498: sampling   200 males and   200 females
Generation   499: sampling   200 males and   200 females
Generation   500: sampling   200 males and   200 females
Generation   501: sampling   200 males and   200 females
In [6]:
ngen,popSize = 1,60
sires1,dams1,gen2 = sampleRan(popSize, ngen, sires1, dams1, gen = gen1);
Generation   502: sampling    30 males and    30 females
In [7]:
ngen,popSize = 1,600
sires1_2,dams1_2,gen2 = sampleRan(popSize, ngen, sires1, dams1, gen = gen1);
Generation   502: sampling   300 males and   300 females
In [8]:
ngen,popSize = 1,1800
siresA1,damsA1,genA1 = sampleRan(popSize, ngen, sires1, dams1_2,gen = gen2)
siresA2,damsA2,genA2 = sampleRan(popSize, ngen, siresA1, damsA1,gen = genA1)
siresA3,damsA3,genA3 = sampleRan(popSize, ngen, siresA2, damsA2,gen = genA2)
siresA4,damsA4,genA4 = sampleRan(popSize, ngen, siresA3, damsA3,gen = genA3)
siresA5,damsA5,genA5 = sampleRan(popSize, ngen, siresA4, damsA4,gen = genA4)
siresA6,damsA6,genA6 = sampleRan(popSize, ngen, siresA5, damsA5,gen = genA5);
Generation   503: sampling   900 males and   900 females
Generation   504: sampling   900 males and   900 females
Generation   505: sampling   900 males and   900 females
Generation   506: sampling   900 males and   900 females
Generation   507: sampling   900 males and   900 females
Generation   508: sampling   900 males and   900 females
In [9]:
animals4 = concatCohorts(siresA4,damsA4)
animals5 = concatCohorts(siresA5,damsA5)
animals6 = concatCohorts(siresA6,damsA6)
#markers
M4 = getOurGenotypes(animals4)
M5 = getOurGenotypes(animals5)
M6 = getOurGenotypes(animals6)
M = [
     M4
     M5
     M6
    ]
#pedigrees
pedigree4 = getPedigree(animals4)
pedigree5 = getPedigree(animals5)
pedigree6 = getPedigree(animals6)
pedigree = [
            pedigree4
            pedigree5
            pedigree6
            ];
In [10]:
k  = size(M,2)
QTLa = 20     # AddQTL
QTLd = 6      # Dominant SNP 
SNPpairs = 3  # SNP–SNP pairs
QTLPos = sample(1:k,QTLa,replace=false)
mrkPos = deleteat!(collect(1:k),sort(QTLPos))
Qa = M[:,sort(QTLPos)]
X = M[:,mrkPos]
newQTLPos = sort(QTLPos) .- (1:QTLa)  # Position of QTL in marker matrix
@printf "The number of additive QTL = %5i  \n" nQTL = size(Qa,2)
@printf "The number of dominant QTL = %5i  \n" QTLd
@printf "The number of epistatic pairs = %5i  \n" SNPpairs
@printf "The maximum MAF = %.2f  \n" maximum(maf(X))
@printf "The number of markers = %5i  \n" nMarkers = size(X,2)
@printf "The number of phenotypic observations = %5i  \n" nObs = size(Qa,1)
The number of additive QTL =    20  
The number of dominant QTL =     6  
The number of epistatic pairs =     3  
The maximum MAF = 0.50  
The number of markers =  3000  
The number of phenotypic observations =  5400  
In [11]:
Xnew = maf(X)
sel1 = vec(Xnew .>= 0.01)
Marker1   = X[:,sel1]
sel2 = vec(Xnew .>= 0.02)
Marker2   = X[:,sel2]
sel3 = vec(Xnew .>= 0.05)
Marker3   = X[:,sel3];
In [12]:
corMat = cor(Marker1)
LDMat = zeros(1800,200)
for i = 1:1800
    LDMat[i,:] = corMat[i,(i+1):(i+200)].^2
end
replace!(LDMat, NaN => 0)
y_Axis = mean(LDMat,dims=1)
plot(x=0.05:0.05:1.0,y=y_Axis[1:20],Guide.title("LD decreases with distance in the first 1cM, MAF>=0.01"), Guide.xlabel("Map Distance,cM"), Guide.ylabel("LD"),Theme(default_color=colorant"blue"))
Out[12]:
Map Distance,cM -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 -1 0 1 2 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 -0.30 -0.28 -0.26 -0.24 -0.22 -0.20 -0.18 -0.16 -0.14 -0.12 -0.10 -0.08 -0.06 -0.04 -0.02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.60 -0.5 0.0 0.5 1.0 -0.30 -0.29 -0.28 -0.27 -0.26 -0.25 -0.24 -0.23 -0.22 -0.21 -0.20 -0.19 -0.18 -0.17 -0.16 -0.15 -0.14 -0.13 -0.12 -0.11 -0.10 -0.09 -0.08 -0.07 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.30 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.40 0.41 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.50 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59 0.60 LD LD decreases with distance in the first 1cM, MAF>=0.01
In [13]:
corMat = cor(Marker2)
LDMat = zeros(1800,200)
for i = 1:1800
    LDMat[i,:] = corMat[i,(i+1):(i+200)].^2
end
replace!(LDMat, NaN => 0)
y_Axis = mean(LDMat,dims=1)
plot(x=0.05:0.05:1.0,y=y_Axis[1:20],Guide.title("LD decreases with distance in the first 1cM, MAF>=0.02"), Guide.xlabel("Map Distance,cM"), Guide.ylabel("LD"),Theme(default_color=colorant"blue"))
Out[13]:
Map Distance,cM -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 -1 0 1 2 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 -0.30 -0.28 -0.26 -0.24 -0.22 -0.20 -0.18 -0.16 -0.14 -0.12 -0.10 -0.08 -0.06 -0.04 -0.02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.60 -0.5 0.0 0.5 1.0 -0.30 -0.29 -0.28 -0.27 -0.26 -0.25 -0.24 -0.23 -0.22 -0.21 -0.20 -0.19 -0.18 -0.17 -0.16 -0.15 -0.14 -0.13 -0.12 -0.11 -0.10 -0.09 -0.08 -0.07 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.30 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.40 0.41 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.50 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59 0.60 LD LD decreases with distance in the first 1cM, MAF>=0.02
In [14]:
corMat = cor(Marker3)
LDMat = zeros(1800,200)
for i = 1:1800
    LDMat[i,:] = corMat[i,(i+1):(i+200)].^2
end
replace!(LDMat, NaN => 0)
y_Axis = mean(LDMat,dims=1)
plot(x=0.05:0.05:1.0,y=y_Axis[1:20],Guide.title("LD decreases with distance in the first 1cM, MAF>=0.05"), Guide.xlabel("Map Distance,cM"), Guide.ylabel("LD"),Theme(default_color=colorant"blue"))
Out[14]:
Map Distance,cM -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 -1 0 1 2 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 -0.30 -0.28 -0.26 -0.24 -0.22 -0.20 -0.18 -0.16 -0.14 -0.12 -0.10 -0.08 -0.06 -0.04 -0.02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.60 -0.5 0.0 0.5 1.0 -0.30 -0.29 -0.28 -0.27 -0.26 -0.25 -0.24 -0.23 -0.22 -0.21 -0.20 -0.19 -0.18 -0.17 -0.16 -0.15 -0.14 -0.13 -0.12 -0.11 -0.10 -0.09 -0.08 -0.07 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.30 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.40 0.41 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.50 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59 0.60 LD LD decreases with distance in the first 1cM, MAF>=0.05
In [15]:
freq = mean(X,dims=1)/2
plot(x=freq,Geom.histogram,
Guide.title("Distribution of Gene Frequency"),
Guide.ylabel("Frequency"),
Guide.xlabel("Gene Frequency"))
Out[15]:
Gene Frequency -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 -1 0 1 2 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 -400 -350 -300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 -500 0 500 1000 -400 -380 -360 -340 -320 -300 -280 -260 -240 -220 -200 -180 -160 -140 -120 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 340 360 380 400 420 440 460 480 500 520 540 560 580 600 620 640 660 680 700 720 740 760 780 800 Frequency Distribution of Gene Frequency
In [16]:
meanXCols = mean(X,dims = 1)
p = meanXCols./2
@printf "The average frequency of heterozygote = %.2f  \n" mean2pq = mean(2p.*(1 .- p))
@printf "The average frequency of p = %.2f  \n" mean_p = mean(p.^2)
@printf "The average frequency of q = %.2f  \n" mean_q = mean((1 .- p).^2)
@printf "Sum of three genotype frequencies = %.2f  \n" mean_p + mean2pq + mean_q
The average frequency of heterozygote = 0.26  
The average frequency of p = 0.37  
The average frequency of q = 0.37  
Sum of three genotype frequencies = 1.00  
In [17]:
sumXCols = sum(X,dims = 1)
geno2sameColnum = size(findall( x -> x == nObs*2, sumXCols))[1]
geno0sameColnum = size(findall( x -> x == 0, sumXCols))[1]
allmumFixLoc = geno2sameColnum + geno0sameColnum
@printf "The number of fixed loci that genotypes are 2 = %5i  \n" geno2sameColnum
@printf "The number of fixed loci that genotypes are 0 = %5i  \n" geno0sameColnum
@printf "Total number of fixed loci = %5i  \n" allmumFixLoc
The number of fixed loci that genotypes are 2 =   236  
The number of fixed loci that genotypes are 0 =   233  
Total number of fixed loci =   469  
In [18]:
α = rand(Gamma(0.4),QTLa) .* sample([-1,1],QTLa);  #scaling parameter α and shape parameter β
In [19]:
d, SortQTLDompos, SortQTLpos = Dominant(QTLPos,QTLd,α);
In [20]:
epistaticQTLPos, epi_coef, epi_combinations = epistatic(QTLPos, SNPpairs);
In [21]:
p = getaf(Qa)
#set fixed allele frequencies to an extreme frequency, so that the algorithm in epi() still works for all loci
p[ p .== 1.0 ] .= 0.99999
p[ p .== 0.0 ] .= 0.00001
#compute genotype frequencies from random mating
gfreq = gf(p,p)
#get epistatic effects
epimodel = "interaction"
epi_vars, epi_effects  = epi(gfreq, p, epi_coef, α, epimodel, epi_combinations);
In [22]:
#get alphas and ds
Newd = zeros(nQTL)
DominQTLpos = findall(x->x in SortQTLDompos, SortQTLpos)
Newd[DominQTLpos] .= d
alphas, ds = get_effects(α,Newd,p,epi_combinations, epi_effects);
Computing orthogonal alphas and ds
In [23]:
#compute noia matrices
noia, noia_dom = make_noia(Qa);
Creating NOIA matrices..
Done creating NOIA matrices!
In [24]:
#breeding values and dominance deviations
simbv = noia * alphas
simdom = noia_dom * ds;
In [25]:
aa = epistaticCom(noia, epi_combinations, noia_dom, nObs, "interaction");
In [26]:
simtgv = simbv + simdom + aa
#phenotypes
resVar = var(simtgv)
resid = rand(Normal(0,sqrt(resVar)), nObs)
pheno = simtgv .+ resid;
In [27]:
h2N = var(simbv)/var(pheno)
Out[27]:
0.5029665893829274
In [28]:
h2B = var(simtgv)/var(pheno)
Out[28]:
0.5089043058964029
In [29]:
QTLPosIndex = [SortQTLpos newQTLPos]
SortQTLDompos = TruePosDom(SortQTLDompos,QTLPosIndex)
sortepistaticQTLPos = TruePosEpis(epistaticQTLPos,QTLPosIndex);
In [30]:
M = Float64.(X)
p = getaf(M)
p[ p .== 1.0 ] .= 0.99999
p[ p .== 0.0 ] .= 0.00001
C = findall(p .< 0.5)
p[C] .= 1 .- p[C]
Z, W = Zmatrix(M, p);
In [31]:
Ka = Z * Z'
Kd = W * W'
Ki = Ka .* Ka - (Z .* Z)*(Z .* Z)';
In [32]:
GA  = (Z'*Z + I*0.00001)./(2*sum(p .*(1 .- p)))
GD = (W'*W + I*0.00001)./(4*sum(p.^2 .* (1 .- p).^2));
GAAtem = GA .* GA + I*0.00001
n = size(GAAtem)[2]
GAA = GAAtem ./ (tr(GAAtem) / n);
In [36]:
ids = string.(1:size(M,1))   
model  = build_model("y = intercept",resVar)
add_genotypes(model,M,var(simbv),header=false,rowID=ids,G_is_marker_variance=false);
The marker IDs are set to 1,2,...,#markers
#markers: 3000; #individuals: 5400
In [37]:
phenTrain = DataFrame(ids=1:size(pheno,1), y=pheno)
out=runMCMC(model, phenTrain,                 # tell JWAS to run analysis, for given model and data 
    Pi=0.99,                                  # intial value of π
    estimatePi=true,
    chain_length=60000,                       # length of chain
    printout_frequency=5000,                  # how often to show progress of analysis 
    printout_model_info=true,                 # tell JWAS to show the options used in this analysis
    methods="BayesC",                         # tell JWAS to run a BayesC analysis
    output_samples_frequency=20,              # how often to output sampled quantities
    output_folder ="results",
    output_samples_for_all_parameters = true,    
    outputEBV=true
);
The folder results is created to save results.
Checking phenotypes...
Individual IDs (strings) are provided in the first column of the phenotypic data.
The number of observations with both genotypes and phenotypes used in the analysis is 5400.

The prior for marker effects variance is calculated from the genetic variance and π.
The mean of the prior for the marker effects variance is: 1.61822



A Linear Mixed Model was build using model equations:

y = intercept

Model Information:

Term            C/F          F/R            nLevels
intercept       factor       fixed                1

MCMC Information:

chain_length                                  60000
burnin                                            0
starting_value                                 true
printout_frequency                             5000
output_samples_frequency                         20
constraint                                    false
missing_phenotypes                             true
update_priors_frequency                           0
seed                                          false

Hyper-parameters Information:

residual variances:                          12.860

Genomic Information:

complete genomic data (i.e., non-single-step analysis)

Genomic Category                               geno
Method                                       BayesC
genetic variances (genomic):                 12.710
marker effect variances:                      1.618
π                                              0.99
estimatePi                                     true
estimateScale                                 false

Degree of freedom for hyper-parameters:

residual variances:                           4.000
marker effect variances:                      4.000



The file results/MCMC_samples_residual_variance.txt is created to save MCMC samples for residual_variance.
The file results/MCMC_samples_marker_effects_geno_y.txt is created to save MCMC samples for marker_effects_geno_y.
The file results/MCMC_samples_marker_effects_variances_geno.txt is created to save MCMC samples for marker_effects_variances_geno.
The file results/MCMC_samples_pi_geno.txt is created to save MCMC samples for pi_geno.
The file results/MCMC_samples_y.intercept.txt is created to save MCMC samples for y:intercept.
The file results/MCMC_samples_EBV_y.txt is created to save MCMC samples for EBV_y.
The file results/MCMC_samples_genetic_variance.txt is created to save MCMC samples for genetic_variance.
running MCMC ...   0%|█                                  |  ETA: 0:04:23
The file results/MCMC_samples_heritability.txt is created to save MCMC samples for heritability.
running MCMC ...   8%|███                                |  ETA: 0:04:26
Posterior means at iteration: 5000
Residual variance: 13.109601
running MCMC ...  17%|██████                             |  ETA: 0:03:55
Posterior means at iteration: 10000
Residual variance: 13.107325
running MCMC ...  17%|██████                             |  ETA: 0:03:55

running MCMC ...  25%|█████████                          |  ETA: 0:03:28
Posterior means at iteration: 15000
Residual variance: 13.108311
running MCMC ...  25%|█████████                          |  ETA: 0:03:28

running MCMC ...  33%|████████████                       |  ETA: 0:03:01
Posterior means at iteration: 20000
Residual variance: 13.109566
running MCMC ...  33%|████████████                       |  ETA: 0:03:01

running MCMC ...  42%|███████████████                    |  ETA: 0:02:38
Posterior means at iteration: 25000
Residual variance: 13.112109
running MCMC ...  42%|███████████████                    |  ETA: 0:02:38

running MCMC ...  50%|██████████████████                 |  ETA: 0:02:16
Posterior means at iteration: 30000
Residual variance: 13.11568
running MCMC ...  50%|██████████████████                 |  ETA: 0:02:16

running MCMC ...  58%|█████████████████████              |  ETA: 0:01:53
Posterior means at iteration: 35000
Residual variance: 13.121654
running MCMC ...  58%|█████████████████████              |  ETA: 0:01:53

running MCMC ...  67%|████████████████████████           |  ETA: 0:01:30
Posterior means at iteration: 40000
Residual variance: 13.12146
running MCMC ...  67%|████████████████████████           |  ETA: 0:01:30

running MCMC ...  75%|███████████████████████████        |  ETA: 0:01:08
Posterior means at iteration: 45000
Residual variance: 13.120212
running MCMC ...  83%|██████████████████████████████     |  ETA: 0:00:45
Posterior means at iteration: 50000
Residual variance: 13.117911
running MCMC ...  92%|█████████████████████████████████  |  ETA: 0:00:23
Posterior means at iteration: 55000
Residual variance: 13.116843
running MCMC ... 100%|███████████████████████████████████|  ETA: 0:00:00
Posterior means at iteration: 60000
Residual variance: 13.115795
running MCMC ... 100%|███████████████████████████████████| Time: 0:04:30

The version of Julia and Platform in use:

Julia Version 1.8.3
Commit 0434deb161 (2022-11-14 20:14 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 4 × Intel(R) Core(TM) i7-7500U CPU @ 2.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, skylake)
  Threads: 1 on 4 virtual cores


The analysis has finished. Results are saved in the returned variable and text files. MCMC samples are saved in text files.


In [38]:
res = GWAS("results/MCMC_samples_marker_effects_geno_y.txt";header=true);
Compute the model frequency for each marker (the probability the marker is included in the model).
In [39]:
maximum(newQTLPos)
Out[39]:
2792
In [40]:
res[newQTLPos,:]
Out[40]:

20 rows × 2 columns

marker_IDmodelfrequency
Abstrac…Float64
12630.003
22850.000333333
33700.00366667
43710.00833333
54240.0
64490.0
76560.00166667
86880.0
97510.0
109900.00366667
1110260.0
1211350.000333333
1315250.424333
1421170.0
1522481.0
1622920.0
1724800.00233333
1826030.0
1926780.000333333
2027920.00433333
In [41]:
out["EBV_y"]
Out[41]:

5,400 rows × 3 columns

IDEBVPEV
AnyAnyAny
11-3.751680.0668002
221.080180.0459823
33-3.948680.0943946
44-0.3394520.0350379
55-5.470860.0522717
660.8233860.0587298
77-4.196250.126959
880.7556640.0350884
99-4.163740.0548528
10100.7847020.125299
1111-3.612960.104193
12124.123920.0440072
1313-0.7812030.0542116
14140.4546960.0342249
1515-4.400290.0697592
1616-0.3892480.0649443
17170.4827670.0488711
18186.038780.0732939
19190.7771650.0355207
2020-1.351450.0678147
2121-6.36020.0535102
22220.1627390.0449542
2323-3.883930.0483218
2424-4.079930.100401
25250.4692090.0329794
2626-0.7897760.0766673
2727-5.49870.0405693
2828-5.481180.0987457
2929-4.362240.0652698
3030-4.237260.0372082
⋮⋮⋮⋮
In [42]:
winVar = GWAS("results/MCMC_samples_marker_effects_geno_y.txt",M;header=true,window_size=150,threshold=0.001)
This function is deprecated. Please check ?GWAS
Out[42]:

20 rows × 6 columns

wStartwEndwSizeprGenVarWPPAPPA_t
Int64Int64Int64Float64Float64Float64
13014501504.011.01.0
221012250150101.881.01.0
3225124001504.521.01.0
4195121001500.770.9996670.999917
5150116501500.440.9316670.986267
6255127001500.340.8026670.955667
7270128501500.050.1650.842714
8135115001500.020.06866670.745958
9285130001500.020.06366670.670148
1090110501500.020.06066670.6092
11240125501500.020.04766670.558152
12165118001500.010.03733330.51475
136017501500.010.03333330.477718
14120113501500.010.02233330.44519
154516001500.00.01633330.4166
161513001500.00.01233330.391333
1711501500.00.0070.368725
187519001500.00.0070.34863
19180119501500.00.005333330.330561
20105112001500.00.003666670.314217
In [43]:
srtIndx = sortperm(winVar.wStart,rev=false)
winVar = DataFrame(wStart = winVar.wStart[srtIndx],
                wEnd   = winVar.wEnd[srtIndx],
                wSize  = winVar.wSize[srtIndx],
                prGenVar = winVar.prGenVar[srtIndx],
                WPPA     = winVar.WPPA[srtIndx],
                PPA_t  = winVar.PPA_t[srtIndx]
            )
Out[43]:

20 rows × 6 columns

wStartwEndwSizeprGenVarWPPAPPA_t
Int64Int64Int64Float64Float64Float64
111501500.00.0070.368725
21513001500.00.01233330.391333
33014501504.011.01.0
44516001500.00.01633330.4166
56017501500.010.03333330.477718
67519001500.00.0070.34863
790110501500.020.06066670.6092
8105112001500.00.003666670.314217
9120113501500.010.02233330.44519
10135115001500.020.06866670.745958
11150116501500.440.9316670.986267
12165118001500.010.03733330.51475
13180119501500.00.005333330.330561
14195121001500.770.9996670.999917
1521012250150101.881.01.0
16225124001504.521.01.0
17240125501500.020.04766670.558152
18255127001500.340.8026670.955667
19270128501500.050.1650.842714
20285130001500.020.06366670.670148
In [44]:
nWins = size(winVar,1)
some_plot1 =plot(x=(1:nWins),y=winVar.WPPA,Guide.title("WPPA of genome windows"), Guide.xlabel("Genomic windows"), Guide.ylabel("WPPA"),
Theme(default_color=colorant"blue"))
Out[44]:
Genomic windows -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 -20 0 20 40 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 -1 0 1 2 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 WPPA WPPA of genome windows
In [45]:
markerEffect=out["marker effects geno"]
add = sortslices([newQTLPos/((numLoci*numChr)-nQTL)*nWins abs.(alphas)],dims=1,by=x->x[1],rev=false)
some_plot2 =plot(layer(x=(1:size(markerEffect,1))*nWins/((numLoci*numChr)-nQTL),y=broadcast(abs,markerEffect.Estimate), Geom.point, Theme(default_color=colorant"blue")),Guide.title("Additive effect SNP scan results"), Guide.xlabel("Genomic windows"), Guide.ylabel("Marker effect"),
layer(x=add[:,1], y=broadcast(abs, add[:,2]), Geom.point, Theme(default_color=colorant"red")))
Out[45]:
Genomic windows -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 -26 -24 -22 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 -25 0 25 50 -25 -24 -23 -22 -21 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 -5 0 5 10 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8 9.0 9.2 9.4 9.6 9.8 10.0 Marker effect Additive effect SNP scan results
In [46]:
DominantLoc = SortQTLDompos
dominant = sortslices([DominantLoc/((numLoci*numChr)-nQTL)*nWins abs.(ds[DominQTLpos])],dims=1,by=x->x[1],rev=false)
some_plot3 =plot(layer(x=(1:size(markerEffect,1))*nWins/((numLoci*numChr)-nQTL),y=broadcast(abs,markerEffect.Estimate), Geom.point, Theme(default_color=colorant"blue")),Guide.title("Dominant effect SNP scan results"), Guide.xlabel("Genomic windows"), Guide.ylabel("Marker effect"),
layer(x=dominant[:,1], y=broadcast(abs, dominant[:,2]), Geom.point, Theme(default_color=colorant"red")))
Out[46]:
Genomic windows -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 -26 -24 -22 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 -25 0 25 50 -25 -24 -23 -22 -21 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 -5 0 5 10 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8 9.0 9.2 9.4 9.6 9.8 10.0 Marker effect Dominant effect SNP scan results
In [47]:
episLocmean = TruePos(epi_combinations,QTLPosIndex)
epistasis = sortslices([episLocmean/((numLoci*numChr)-nQTL)*nWins abs.(sum(epi_effects[:,[5]],dims=2))],dims=1,by=x->x[1],rev=false)
some_plot4 =plot(layer(x=(1:size(markerEffect,1))*nWins/((numLoci*numChr)-nQTL),y=broadcast(abs,markerEffect.Estimate), Geom.point, Theme(default_color=colorant"blue")),Guide.title("Epistatic effect SNP scan results"), Guide.xlabel("Genomic windows"), Guide.ylabel("Marker effect"),
layer(x=epistasis[:,1], y=broadcast(abs, epistasis[:,2]), Geom.point, Theme(default_color=colorant"red")))
Out[47]:
Genomic windows -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 -26 -24 -22 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 -25 0 25 50 -25 -24 -23 -22 -21 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 -5 0 5 10 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8 9.0 9.2 9.4 9.6 9.8 10.0 Marker effect Epistatic effect SNP scan results
In [48]:
out["marker effects geno"]
Out[48]:

3,000 rows × 5 columns

TraitMarker_IDEstimateSDModel_Frequency
AnyAnyAnyAnyAny
1y10.00.00.0
2y21.28741f-50.0007050230.000333334
3y30.00.00.0
4y40.00.00.0
5y50.00.00.0
6y60.00.00.0
7y7-0.005009980.1547230.00366667
8y80.00.00.0
9y9-4.83103f-50.001880280.000666667
10y100.0004437060.1384070.00366667
11y11-0.0001159060.003970870.001
12y120.00.00.0
13y130.00.00.0
14y140.0001141460.00456430.000666668
15y150.00.00.0
16y160.003780670.05933660.00466667
17y170.00.00.0
18y18-7.0054f-50.003836380.000333333
19y190.00.00.0
20y200.002714040.1495640.00533334
21y210.00.00.0
22y220.0001673740.005713940.001
23y230.0001976190.006442860.001
24y240.0003079760.009860870.001
25y250.00.00.0
26y260.0001650460.1479040.00333333
27y270.00.00.0
28y280.00.00.0
29y293.1048f-50.001700280.000333334
30y300.002026140.1219840.00266667
⋮⋮⋮⋮⋮⋮
In [49]:
wSize = 150
wStartV,wEndV,testStat = winTest(M[1:500,:],pheno[1:500],var(simbv),resVar,wSize)
srtIndx = sortperm(testStat,rev=true)
replace!(testStat, NaN => 0)
bigF = [wStartV[srtIndx] wEndV[srtIndx] testStat[srtIndx]];
some_plot5 =plot(x=(1:nWins),y=testStat,Guide.title("F statistics of Genome window "), Guide.xlabel("Genomic windows"), Guide.ylabel("F value"),
Theme(default_color=colorant"blue"))
Out[49]:
Genomic windows -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 -20 0 20 40 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 -5 0 5 10 15 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8 9.0 9.2 9.4 9.6 9.8 10.0 10.2 10.4 10.6 10.8 11.0 F value F statistics of Genome window
In [50]:
sortPosQTL = sort(QTLPos)
wPos = [bigF [findfirst(bigF[i,1] .<= sortPosQTL .< bigF[i,2]) for i=1:size(bigF,1)]]
Out[50]:
20×4 Matrix{Union{Nothing, Float64}}:
 1051.0  1200.0  6.50106  12.0
 2551.0  2700.0  5.35068  18.0
 2851.0  3000.0  5.23125    nothing
 1801.0  1950.0  5.18532    nothing
 2701.0  2850.0  4.86845  20.0
 1951.0  2100.0  4.73664    nothing
  301.0   450.0  4.70326   3.0
 1201.0  1350.0  4.63476    nothing
 1501.0  1650.0  4.63083  13.0
  601.0   750.0  4.62811   7.0
  151.0   300.0  4.62026   1.0
  901.0  1050.0  4.2378   10.0
 2101.0  2250.0  4.17256  14.0
    1.0   150.0  4.04736    nothing
 1351.0  1500.0  3.90571    nothing
  751.0   900.0  3.87372   9.0
 2401.0  2550.0  3.78504  17.0
 2251.0  2400.0  3.69143  15.0
  451.0   600.0  3.52503   6.0
 1651.0  1800.0  3.2208     nothing
In [ ]: